/*This script generates results at the PLSS level data found in the main text;
Table 4
Figure 8
Figure 9

Table D1 (tabular support for Figure 8)
Table D7 (tabular support for Figure 9)


*/


cd "$data"
set more off

*Bring in the PLSS data
use PLSS_Turbines_Irrigation_JAERE.dta, replace

*Globals drawn on
global xcontrols elev_mean elev_std soil stream3_dist tmean ppt_mean lat lon 
global xend sharecropland sharedeveloped shareforest powerlines_dist
global xwind ave_wind_class max_wind_class
global SE vce(cl alt_county_id)  


/*Table 4*/
local i = 1
foreach irr in "irr_mean" "CPIS_per" "irr_mean CPIS_per" {
	 reghdfe turbine `irr' $xwind $xcontrols $xend , absorb(alt_county_id) $SE
		qui sum turbine if e(sample) ==1 
		local mdv = r(mean)
		estadd scalar MDV = `mdv'
		estadd scalar Nfe = e(df_a_initial)
		estadd scalar Nclus =  e(N_clust1)
	
	estimates store turb_irr_cpis_`i'

	local i = `i' + 1
	}

foreach irr in "irr_mean" "CPIS_per"  {
	qui logit turbine `irr' $xwind $xcontrols $xend i.alt_county_id if turb_alt_county_id_ind==1, $SE
		qui sum turbine if e(sample) ==1 
		local mdv = r(mean)
		margins, dydx(`irr') 
		local odd_`irr'=el(r(b),1,1)
		estadd scalar MDV = `mdv'
		estadd scalar Margins_`irr'= `odd_`irr''
		estadd scalar Nfe = e(k)-16
		estadd scalar Nclus =  e(N_clust)
	
	estimates store turb_irr_cpis_`i'
	
	
	local i=`i'+1		
	}
	
qui logit turbine irr_mean CPIS_per $xwind $xcontrols $xend i.alt_county_id if turb_alt_county_id_ind==1, $SE
	qui sum turbine if e(sample) ==1 
	local mdv = r(mean)
	margins, dydx(irr_mean) 
	local odd_irr_mean=el(r(b),1,1)
	margins, dydx(CPIS_per) 
	local odd_CPIS_per=el(r(b),1,1)
	estadd scalar MDV = `mdv'
	estadd scalar Margins_irr_mean = `odd_irr_mean'
	estadd scalar Margins_CPIS_per = `odd_CPIS_per'
	estadd scalar Nfe = e(k)-17
	estadd scalar Nclus =  e(N_clust)
	
estimates store turb_irr_cpis_`i'

esttab turb_irr_cpis_* using "$results/tables/tab4_PLSS_turb_CPIS_irr.csv",  label  keep(irr_mean irrigated CPIS_per cpis_ind2 ave_wind_class max_wind_class sharecropland) order(irr_mean irrigated CPIS_per cpis_ind2 sharecropland ave_wind_class max_wind_class) star(* 0.1 ** 0.05 *** 0.01)  ar2 pr2 se(a3) b(a3) scalars("MDV Mean DV" "Margins_irr_mean MEM of Irrigation" "Margins_CPIS_per MEM of CPIS" "Nfe N Fixed Effects" "Nclus Standard Error Clusters") replace


/*Figure 8 and Table D1*/
		
local s = 0	
foreach irr in "irr_mean" "CPIS_per" "irr_mean CPIS_per" {	
	local s = `s' + 1
	local i = 1
 	foreach fe in alt_state_id alt_county_id alt_plssid fe_town_3 fe_town_4 fe_town_9 {
		
		qui reghdfe turbine `irr' $xwind $xcontrols $xend , absorb(`fe') $SE
			qui sum turbine if e(sample) ==1 
			local mdv = r(mean)
			estadd scalar MDV = `mdv'
			estadd scalar Nfe = e(df_a_initial)
			estadd scalar Nclus =  e(N_clust1)
	
		estimates store turb_fe_rob_`s'_`i'

		local i=`i'+1
	
		}
	}	

*Each panel (A,B, and C) is exported separately*
	
esttab turb_fe_rob_1_* using "$results/tables/tabD1a_PLSS_turb_FE.csv",  label  keep(irr_mean CPIS_per ) order(irr_mean CPIS_per ) star(* 0.1 ** 0.05 *** 0.01)  ar2 se(a3) b(a3) scalars("MDV Mean DV"  "Nfe N Fixed Effects" "Nclus Standard Error Clusters") replace

esttab turb_fe_rob_2_* using "$results/tables/tabD1b_PLSS_turb_FE.csv",  label  keep(irr_mean CPIS_per ) order(irr_mean CPIS_per ) star(* 0.1 ** 0.05 *** 0.01)  ar2 se(a3) b(a3) scalars("MDV Mean DV"  "Nfe N Fixed Effects" "Nclus Standard Error Clusters") replace

esttab turb_fe_rob_3_* using "$results/tables/tabD1c_PLSS_turb_FE.csv",  label  keep(irr_mean CPIS_per ) order(irr_mean CPIS_per ) star(* 0.1 ** 0.05 *** 0.01)  ar2 se(a3) b(a3) scalars("MDV Mean DV"  "Nfe N Fixed Effects" "Nclus Standard Error Clusters") replace

*plot of Fixed Effect Specifications
preserve
	
		
		clear
		local binlist  1 2 3 4 5 6
		local numbins = wordcount("`binlist'")
			
		set obs `numbins'
		gen bin =.
			
		
		foreach access of numlist 0 1  {
			gen beta_type`access' =.
			gen se_type`access' =.	
			}

			gen id = _n
				
				
			forv i = 1/`numbins' {
				local bin: word `i' of `binlist'
				replace bin = `bin' if id == `i'
				
				estimates restore  turb_fe_rob_3_`bin'

						
				replace beta_type0 = _b[irr_] if id == `i'
				replace se_type0 = _se[irr_] if id == `i'
		
				lincom _b[irr_] +_b[CPIS_per]

				replace beta_type1 = r(estimate) if id == `i'
				replace se_type1 = r(se) if id == `i'
			
		
				}
		
				gen plus90_0=beta_type0+1.64*se_type0
				gen minus90_0=beta_type0-1.64*se_type0
				gen plus90_1=beta_type1+1.64*se_type1
				gen minus90_1=beta_type1-1.64*se_type1
		
			
				replace bin=bin-.05
				gen bin_alt=bin+0.1
				
		label define fixed_effect 1 "State" 2 "County" 3 "Township" 4 "Township Thirds" 5 "Township Fourths" 6 "Township Ninths"
		label val bin fixed_effect

		label var bin "Spatial Fixed Effect Level"
		
	
		
		twoway (scatter beta_type0 bin, lc(blue) mc(blue) m(S)) (rcap plus90_0 minus90_0 bin, lpattern(dash) lcolor(blue) msize(large)) ///
		(scatter beta_type1 bin_alt, lc(forest_green) mc(forest_green) m(O)) (rcap plus90_1 minus90_1 bin_alt, lpattern(solid) lcolor(forest_green) msize(large)), ///
		 graphregion(fcolor(white) style(none) color(gs16)) ytitle("{&Delta}Pr(Turbine)" ) xtitle("Spatial Fixed Effect Level") xlabel(,valuelabel angle(-45)) xscale(range(0.5 6.5)) ylabel( , ) yline(0) legend(label(1 "Non-CPIS Irrigation") label(3 "CPIS Irrigation") order(1 3) region(lstyle(none))) note("90th% Confidence Interval") xsize(6.5) ysize(3.5)
				graph export  "$results/figures/fig8_plss_CPISIrr_FE_bin.pdf", replace
		
		restore

/*Figure 9 and Table D7*/

local i = 1
foreach combo in "CPIS_per othercpis_dist" "CPIS_per cpis_neighbors" "CPIS_per i.cpis_neighbors" "i.cpis_ind2 i.cpis_neighbors" "i.CPIS_bin i.cpis_neighbors" "cpis_ind2##i.cpis_neighbors" {
 
	qui reghdfe turbine `combo' irrigated $xwind $xcontrols $xend , absorb(alt_county_id) $SE
		qui sum turbine if e(sample) ==1 
		local mdv = r(mean)
		estadd scalar MDV = `mdv'
		estadd scalar Nfe = e(df_a_initial)
		estadd scalar Nclus =  e(N_clust1)
	
	estimates store turb_neighbor_`i'
	
	local i = `i'+1
	}

foreach rest in 0 1 {
	reghdfe turbine irrigated i.cpis_neighbors $xwind $xcontrols $xend if cpis_ind2==`rest', absorb(alt_county_id) $SE
 	
		qui sum turbine if e(sample) ==1 
		local mdv = r(mean)
		estadd scalar MDV = `mdv'
		estadd scalar Nfe = e(df_a_initial)
		estadd scalar Nclus =  e(N_clust1)
	
	estimates store turb_neighbor_`i'
	
	local i = `i'+1
	}
	
*Column (8) is rearranged into two columns in the print version
esttab turb_neighbor_1 turb_neighbor_2 turb_neighbor_3 turb_neighbor_4 turb_neighbor_5 turb_neighbor_7 turb_neighbor_8 turb_neighbor_6  using "$results/tables/tabD7_PLSS_turb_CPIS_neighbor_irr.csv",  label  keep(1.CPIS_bin 2.CPIS_bin 3.CPIS_bin 4.CPIS_bin 5.CPIS_bin CPIS_per 1.cpis_ind2 irrigated othercpis_dist cpis_neighbors 1.cpis_neighbors 2.cpis_neighbors 3.cpis_neighbors 4.cpis_neighbors 5.cpis_neighbors 6.cpis_neighbors 7.cpis_neighbors 8.cpis_neighbors 1.cpis_ind2#1.cpis_neighbors 1.cpis_ind2#2.cpis_neighbors 1.cpis_ind2#3.cpis_neighbors 1.cpis_ind2#4.cpis_neighbors 1.cpis_ind2#5.cpis_neighbors 1.cpis_ind2#6.cpis_neighbors 1.cpis_ind2#7.cpis_neighbors 1.cpis_ind2#8.cpis_neighbors) order(othercpis_dist cpis_neighbors 1.cpis_neighbors 2.cpis_neighbors 3.cpis_neighbors 4.cpis_neighbors 5.cpis_neighbors 6.cpis_neighbors 7.cpis_neighbors 8.cpis_neighbors irrigated CPIS_per 1.cpis_ind2 1.CPIS_bin 2.CPIS_bin 3.CPIS_bin 4.CPIS_bin 5.CPIS_bin ) star(* 0.1 ** 0.05 *** 0.01)  ar2 se(a3) b(a3) scalars("MDV Mean DV" "Nfe N Fixed Effects" "Nclus Standard Error Clusters") replace


**plot of neighbor effects
preserve
	
		
		clear
		estimates restore  turb_neighbor_6
		local binlist  0 1 2 3 4 5 6 7 8
		local numbins = wordcount("`binlist'")
			
		set obs `numbins'
		gen bin =.
		
		foreach c in 0 1 {
			gen beta`c'=.
			gen se`c'=.
			}
		
		gen id=_n
		
		forvalue i = 1/`numbins' {
			local bin: word `i' of `binlist'
			replace bin = `bin' if id == `i'
			
			replace beta0=_b[`bin'.cpis_neighbors] if id == `i'
			
			replace se0=_se[`bin'.cpis_neighbors] if id == `i'
			
			lincom [1.cpis_ind2#`bin'.cpis_neighbors+1.cpis_ind2+`bin'.cpis_neighbor] 
			
			replace beta1 = r(estimate) if id == `i'
			replace se1 = r(se) if id==`i'
			
			}
		
		
		foreach c in 0 1 {
			gen plus95_`c'=beta`c'+1.96*se`c'
			gen minus95_`c'=beta`c'-1.96*se`c'

			gen plus90_`c'=beta`c'+1.64*se`c'
			gen minus90_`c'=beta`c'-1.64*se`c'
			}
			


		gen bin0=bin-.1
		gen bin1=bin+.1
	
			label var bin1 "Neighboring Sections with CPIS"
			label var bin0 "Neighboring Sections with CPIS"
			label var bin "Neighboring Sections with CPIS"

		
		twoway (scatter beta0 bin0, lc(maroon) mc(maroon) msymbol(D)) (rcap plus90_0 minus90_0 bin0, lpattern(dash) lcolor(maroon) msize(large))  ///
		 (scatter beta1 bin1, lc(forest_green) mc(forest_green) msymbol(C)) (rcap plus90_1 minus90_1 bin1, lpattern(solid) lcolor(forest_green) msize(large)) , ///
		scheme(sj) graphregion(fcolor(white) style(none) color(gs16)) ytitle("{&Delta}Pr(Turbine)" )  xlabel(0(1)8) ylabel( , nogrid) yline(0) legend(label(1 "Sections w/o CPIS") label(3 "Sections w/ CPIS") order(1 3) region(lstyle(none))) note("90th% Confidence Interval") xsize(6.5) ysize(3.5)
				graph export  "$results/figures/fig9_plss_CPIS_neighbors_bif_irr.pdf", replace
		
		restore 

****This completes the PLSS analysis in the main text
